! ******************************************************
! Boundary Conditions
! ******************************************************

module Boundary
    use Initialization
    use Injector
    implicit none

contains

    subroutine applyParticleBoundary()
        implicit none
        if(UseInjector) call applyInjectorForParticle()
        call applyOpenForParticle()
        call applyPeriodicForParticle()
        call applyReflectForParticle()
        call applyThermalForParticle()
        
    end subroutine applyParticleBoundary

    subroutine applyGridBoundary()
        implicit none

        call applyOpenForGrid()
        call applyReflectForGrid()
        call applyPeriodicForGrid()
        call applyThermalForGrid()
    end subroutine applyGridBoundary

    subroutine applyGridBoundary_in_exchange(G)
        implicit none
        real(kind=8) :: G(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)

        call applyOpenForGrid_in_exchange(G)
        call applyReflectForGrid_in_exchange(G)
        call applyPeriodicForGrid_in_exchange(G)
        call applyThermalForGrid_in_exchange(G)
    end subroutine applyGridBoundary_in_exchange

    subroutine applyOpenForGrid()
        implicit none
        if(UseInjector)then
            call usePeriodicForGrid_Injector_for(ni)
            call usePeriodicForGrid_Injector_for(ne)
            call usePeriodicForGrid_Injector_for(ux)
            call usePeriodicForGrid_Injector_for(uy)
            call usePeriodicForGrid_Injector_for(uz)
        else
            call useOpen3ForGrid(1)
        end if
    end subroutine applyOpenForGrid

    subroutine applyOpenForGrid_in_exchange(G)
        implicit none
        real(kind=8) :: G(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)
        if(UseInjector)then
            call usePeriodicForGrid_Injector_in_exchange(G)
        else
            call useOpen3ForGrid_in_exchange(1,G)   
        end if
    end subroutine applyOpenForGrid_in_exchange

    subroutine applyPeriodicForGrid()
        implicit none
        call applyPeriodicForGridfor(ni)
        call applyPeriodicForGridfor(ne)
        call applyPeriodicForGridfor(ux)
        call applyPeriodicForGridfor(uy)
        call applyPeriodicForGridfor(uz)
        ! call applyPeriodicForGridfor(Ex)
        ! call applyPeriodicForGridfor(Ey)
        ! call applyPeriodicForGridfor(Ez)
        ! call applyPeriodicForGridfor(Bx)
        ! call applyPeriodicForGridfor(By)
        ! call applyPeriodicForGridfor(Bz)

    end subroutine applyPeriodicForGrid

    subroutine applyPeriodicForGrid_in_exchange(G)
        implicit none
        real(kind=8) :: G(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)
        
        call applyPeriodicForGrid_in_exchange_for(G)

    end subroutine applyPeriodicForGrid_in_exchange

    subroutine applyReflectForGrid()
        implicit none

        
        call useReflectForGrid(3)
        ! call useXinverseForGrid(3)
    end subroutine applyReflectForGrid

    subroutine applyReflectForGrid_in_exchange(G)
        implicit none
        real(kind=8) :: G(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)

        
        call useReflectForGrid_in_exchange(3,G)
        ! call useXinverseForGrid_in_exchange(3,G)
    end subroutine applyReflectForGrid_in_exchange

    subroutine applyThermalForGrid()
        implicit none

        call useOpen3ForGrid(4)
    end subroutine applyThermalForGrid

    subroutine applyThermalForGrid_in_exchange(G)
        implicit none
        real(kind=8) :: G(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)

        call useOpen3ForGrid_in_exchange(4,G)
    end subroutine applyThermalForGrid_in_exchange

    ! Open Condition
    subroutine applyOpenForParticle()
        implicit none

        if(xMinBoundary==1)then
            if(ppi==0)then
                call deleteList(MPIListHeadLeft)
            end if
        end if

        if(xMaxBoundary==1)then
            if(ppi==ppx-1)then
                call deleteList(MPIListHeadRight)
            end if
        end if

        if(yMinBoundary==1)then
            if(ppj==0)then
                call deleteList(MPIListHeadDown)
            end if
        end if

        if(yMaxBoundary==1)then
            if(ppj==ppy-1)then
                call deleteList(MPIListHeadUp)
            end if
        end if
    end subroutine applyOpenForParticle

    ! Periodic Condition
    subroutine applyPeriodicForParticle()
        implicit none
        integer(kind=4) :: obj

        if(xMinBoundary==2 .or. xMaxBoundary==2)then
            if(ppi==0)then
                obj = ppx*(ppj+1)-1
                call setPeriodicPosition(MPIListHeadLeft)
                call sendIons(obj,MPIListHeadLeft)
                call recvIons(obj,IonListHead)
            else if(ppi==ppx-1)then
                obj = ppx*ppj
                call recvIons(obj,IonListHead)
                call setPeriodicPosition(MPIListHeadRight)
                call sendIons(obj,MPIListHeadRight)
            end if
        end if

        if(yMinBoundary==2 .or. yMaxBoundary==2)then
            if(ppj==0)then
                obj = ppx*(ppy-1)+ppi
                call setPeriodicPosition(MPIListHeadDown)
                call sendIons(obj,MPIListHeadDown)
                call recvIons(obj,IonListHead)
            else if(ppj==ppy-1)then
                obj = ppi
                call recvIons(obj,IonListHead)
                call setPeriodicPosition(MPIListHeadUp)
                call sendIons(obj,MPIListHeadUp)
            end if
        end if
    end subroutine applyPeriodicForParticle

    subroutine applyPeriodicForGridfor(G)
        implicit none
        integer(kind=4) :: obj
        integer(kind=4) :: num, size
        real(kind=8) :: Tempx(2+ncx,-1-ncy:Ny+ncy+1), Tempy(-1-ncx:Nx+ncx+1,2+ncy)
        real(kind=8) :: G(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)
        size = 2+ncx
        
        if(xMinBoundary==2 .or. xMaxBoundary==2)then
            num = Ny+3+2*ncy
            if(ppi==0)then
                obj = ppx*(ppj+1)-1
                Tempx(1:2+ncx,:) = G(0:1+ncx,:)
                call MPI_SEND(Tempx,size*num,mpi_real8,obj,1,mpi_comm_world,ierr)
                call MPI_RECV(Tempx,size*num,mpi_real8,obj,1,mpi_comm_world, &
                              status,ierr)
                G(-1-ncx:0,:) = Tempx(1:2+ncx,:)
            else if(ppi==ppx-1)then
                obj = ppx*ppj
                call MPI_RECV(Tempx,size*num,mpi_real8,obj,1,mpi_comm_world, &
                              status,ierr)
                if(yMinBoundary==1 .and. ppj==0)then
                    G(Nx,0:Ny+ncy+1) = G(Nx,0:Ny+ncy+1)+Tempx(1,0:Ny+ncy+1)
                elseif(yMaxBoundary==1 .and. ppj==ppy-1)then
                    G(Nx,-1-ncy:Ny) = G(Nx,-1-ncy:Ny)+Tempx(1,-1-ncy:Ny)
                else
                    G(Nx,:) = G(Nx,:)+Tempx(1,:)
                end if
                
                G(Nx+1:Nx+ncx+1,:) = Tempx(2:2+ncx,:)
                
                Tempx(1:2+ncx,:) = G(Nx-ncx-1:Nx,:)
                call MPI_SEND(Tempx,size*num,mpi_real8,obj,1,mpi_comm_world,ierr)
            end if
        end if
        size = 2+ncy
        if(yMinBoundary==2 .or. yMaxBoundary==2)then
            num = Nx+3+2*ncx
            if(ppj==0)then
                obj = ppx*(ppy-1)+ppi
                Tempy(:,1:2+ncy) = G(:,0:1+ncy)
                call MPI_SEND(Tempy,num*size,mpi_real8,obj,1,mpi_comm_world,ierr)
                call MPI_RECV(Tempy,num*size,mpi_real8,obj,1,mpi_comm_world, &
                              status,ierr)
                G(:,-1-ncy:0) = Tempy(:,1:2+ncy)
            else if(ppj==ppy-1)then
                obj = ppi
                call MPI_RECV(Tempy,num*size,mpi_real8,obj,1,mpi_comm_world, &
                              status,ierr)
                if(xMinBoundary==1 .and. ppi==0)then
                    G(0:Nx+ncx+1,Ny) = G(0:Nx+ncx+1,Ny)+Tempy(0:Nx+ncx+1,1)
                elseif(xMaxBoundary==1 .and. ppi==ppx-1)then
                    G(-1-ncx:Nx,Ny) = G(-1-ncx:Nx,Ny)+Tempy(-1-ncx:Nx,1)
                else
                    G(:,Ny) = G(:,Ny)+Tempy(:,1)
                end if
                G(:,Ny+1:Ny+ncy+1) = Tempy(:,2:2+ncy)

                Tempy(:,1:2+ncy) = G(:,Ny-ncy-1:Ny)
                call MPI_SEND(Tempy,num*size,mpi_real8,obj,1,mpi_comm_world,ierr)
            end if
        end if
    end subroutine applyPeriodicForGridfor

    subroutine applyPeriodicForGrid_in_exchange_for(G)
        implicit none
        integer(kind=4) :: obj
        integer(kind=4) :: num, size
        real(kind=8) :: Tempx(2+ncx,-1-ncy:Ny+ncy+1), Tempy(-1-ncx:Nx+ncx+1,2+ncy)
        real(kind=8) :: Tempx_b(1+ncx,-1-ncy:Ny+ncy+1), Tempy_b(-1-ncx:Nx+ncx+1,1+ncy)
        real(kind=8) :: G(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)
        
        if(xMinBoundary==2 .or. xMaxBoundary==2)then
            if(loc(G)==loc(By) .or. loc(G)==loc(Bz) .or. &
            loc(G)==loc(Bfy) .or. loc(G)==loc(Bfz) .or. &
            loc(G)==loc(Bwy) .or. loc(G)==loc(Bwz) .or. loc(G)==loc(Ex))then
                size = 1+ncx
                num = Ny+3+2*ncy
                if(ppi==0)then
                    obj = ppx*(ppj+1)-1
                    Tempx_b(1:1+ncx,:) = G(1:1+ncx,:)
                    call MPI_SEND(Tempx_b,size*num,mpi_real8,obj,1,mpi_comm_world,ierr)
                    call MPI_RECV(Tempx_b,size*num,mpi_real8,obj,1,mpi_comm_world, &
                                status,ierr)
                    G(-1-ncx:-1,:) = Tempx_b(1:1+ncx,:)
                else if(ppi==ppx-1)then
                    obj = ppx*ppj
                    call MPI_RECV(Tempx_b,size*num,mpi_real8,obj,1,mpi_comm_world, &
                                status,ierr)
                    
                    G(Nx+1:Nx+ncx+1,:) = Tempx_b(1:1+ncx,:)
                    
                    Tempx_b(1:1+ncx,:) = G(Nx-ncx:Nx,:)
                    call MPI_SEND(Tempx_b,size*num,mpi_real8,obj,1,mpi_comm_world,ierr)
                end if
            else
                size = 2+ncx
                num = Ny+3+2*ncy
                if(ppi==0)then
                    obj = ppx*(ppj+1)-1
                    Tempx(1:2+ncx,:) = G(0:1+ncx,:)
                    call MPI_SEND(Tempx,size*num,mpi_real8,obj,1,mpi_comm_world,ierr)
                    call MPI_RECV(Tempx,size*num,mpi_real8,obj,1,mpi_comm_world, &
                                status,ierr)
                    G(-1-ncx:0,:) = Tempx(1:2+ncx,:)
                else if(ppi==ppx-1)then
                    obj = ppx*ppj
                    call MPI_RECV(Tempx,size*num,mpi_real8,obj,1,mpi_comm_world, &
                                status,ierr)
                    G(Nx,:) = (G(Nx,:)+Tempx(1,:))/2
                    G(Nx+1:Nx+ncx+1,:) = Tempx(2:2+ncx,:)
                    
                    Tempx(1:2+ncx,:) = G(Nx-ncx-1:Nx,:)
                    call MPI_SEND(Tempx,size*num,mpi_real8,obj,1,mpi_comm_world,ierr)
                end if
            end if
        end if
        
        if(yMinBoundary==2 .or. yMaxBoundary==2)then
            if(loc(G)==loc(Bx) .or. loc(G)==loc(Bz) .or. &
            loc(G)==loc(Bfx) .or. loc(G)==loc(Bfz) .or. &
            loc(G)==loc(Bwx) .or. loc(G)==loc(Bwz) .or. loc(G)==loc(Ey))then
                size = 1+ncy
                num = Nx+3+2*ncx
                if(ppj==0)then
                    obj = ppx*(ppy-1)+ppi
                    Tempy_b(:,1:1+ncy) = G(:,1:1+ncy)
                    call MPI_SEND(Tempy_b,num*size,mpi_real8,obj,1,mpi_comm_world,ierr)
                    call MPI_RECV(Tempy_b,num*size,mpi_real8,obj,1,mpi_comm_world, &
                                status,ierr)
                    G(:,-1-ncy:-1) = Tempy_b(:,1:1+ncy)
                else if(ppj==ppy-1)then
                    obj = ppi
                    call MPI_RECV(Tempy_b,num*size,mpi_real8,obj,1,mpi_comm_world, &
                                status,ierr)
                    
                    G(:,Ny+1:Ny+ncy+1) = Tempy_b(:,1:1+ncy)

                    Tempy_b(:,1:1+ncy) = G(:,Ny-ncy:Ny)
                    call MPI_SEND(Tempy_b,num*size,mpi_real8,obj,1,mpi_comm_world,ierr)
                end if
            else
                size = 2+ncy
                num = Nx+3+2*ncx
                if(ppj==0)then
                    obj = ppx*(ppy-1)+ppi
                    Tempy(:,1:2+ncy) = G(:,0:1+ncy)
                    call MPI_SEND(Tempy,num*size,mpi_real8,obj,1,mpi_comm_world,ierr)
                    call MPI_RECV(Tempy,num*size,mpi_real8,obj,1,mpi_comm_world, &
                                status,ierr)
                    G(:,-1-ncy:0) = Tempy(:,1:2+ncy)
                else if(ppj==ppy-1)then
                    obj = ppi
                    call MPI_RECV(Tempy,num*size,mpi_real8,obj,1,mpi_comm_world, &
                                status,ierr)
                    
                    G(:,Ny) = (G(:,Ny)+Tempy(:,1))/2
                    G(:,Ny+1:Ny+ncy+1) = Tempy(:,2:2+ncy)

                    Tempy(:,1:2+ncy) = G(:,Ny-ncy-1:Ny)
                    call MPI_SEND(Tempy,num*size,mpi_real8,obj,1,mpi_comm_world,ierr)
                end if
            end if
        end if
    end subroutine applyPeriodicForGrid_in_exchange_for

    ! Reflect Condition
    subroutine applyReflectForParticle()
        implicit none
        type(PseudoParticle), pointer :: p, q
        logical :: IonIsInside

        if(xMinBoundary==3)then
            if(ppi==0)then
                q => MPIListHeadLeft
                p => MPIListHeadLeft%next
                do while(associated(p))
                    p%x = abs(p%x-xMinBox)+xMinBox
                    p%vx = -p%vx
                    call isInside(p,IonIsInside)
                    if(IonIsInside)then
                        call insertIon(IonListHead,p)
                        call deleteIon(MPIListHeadLeft,q)
                        p => q%next
                    else
                        call deleteIon(MPIListHeadLeft,q)
                        p => q%next
                    end if
                end do
            end if
        end if

        if(xMaxBoundary==3)then
            if(ppi==ppx-1)then
                q => MPIListHeadRight
                p => MPIListHeadRight%next
                do while(associated(p))
                    p%x = xMaxBox-abs(p%x-xMaxBox)
                    p%vx = -p%vx
                    call isInside(p,IonIsInside)
                    if(IonIsInside)then
                        call insertIon(IonListHead,p)
                        call deleteIon(MPIListHeadRight,q)
                        p => q%next
                    else
                        call deleteIon(MPIListHeadRight,q)
                        p => q%next
                    end if
                end do
            end if
        end if

        if(yMinBoundary==3)then
            if(ppj==0)then
                q => MPIListHeadDown
                p => MPIListHeadDown%next
                do while(associated(p))
                    p%y = abs(p%y-yMinBox)+yMinBox
                    p%vy = -p%vy
                    call isInside(p,IonIsInside)
                    if(IonIsInside)then
                        call insertIon(IonListHead,p)
                        call deleteIon(MPIListHeadDown,q)
                        p => q%next
                    else
                        call deleteIon(MPIListHeadDown,q)
                        p => q%next
                    end if
                end do
            end if
        end if

        if(yMaxBoundary==3)then
            if(ppj==ppy-1)then
                q => MPIListHeadUp
                p => MPIListHeadUp%next
                do while(associated(p))
                    p%y = yMaxBox-abs(p%y-yMaxBox)
                    p%vy = -p%vy
                    call isInside(p,IonIsInside)
                    if(IonIsInside)then
                        call insertIon(IonListHead,p)
                        call deleteIon(MPIListHeadUp,q)
                        p => q%next
                    else
                        call deleteIon(MPIListHeadUp,q)
                        p => q%next
                    end if
                end do
            end if
        end if
    end subroutine applyReflectForParticle

    ! Thermal Condition
    subroutine applyThermalForParticle()
        implicit none
        type(PseudoParticle), pointer :: EdgeListHead, CornerListHead
        type(PseudoParticle), pointer :: p, q
        integer(kind=4) :: i, j, n
        logical :: CellIsEmpty

        call createHeadOf(EdgeListHead)
        call createHeadOf(CornerListHead)

        ! Edge Ions
        if(xMinBoundary==4 .and. ppi==0)then
            call deleteList(MPIListHeadLeft)
            q => EdgeListHead
            do n=1,NSpecies
                do j=0,Ny-1
                    CellIsEmpty = isCellEmpty(-1,j,n)
                    if(.not. CellIsEmpty)then
                        call createIons(-1,j,n,q,EdgeListHead)
                    end if
                end do
            end do

            q => EdgeListHead
            p => EdgeListHead%next
            do while(associated(p))
                p%x = p%x+p%vx*dt
                p%y = p%y+p%vy*dt
                if(p%x>xMinBox .and. p%y>yMinBox .and. p%y<yMaxBox)then
                    q => p
                    p => p%next
                else
                    call deleteIon(EdgeListHead,q)
                    p => q%next
                end if
            end do
            call moveListA2B(EdgeListHead,IonListHead)
        end if

        if(xMaxBoundary==4 .and. ppi==ppx-1)then
            call deleteList(MPIListHeadRight)
            q => EdgeListHead
            do n=1,NSpecies
                do j=0,Ny-1
                    CellIsEmpty = isCellEmpty(Nx,j,n)
                    if(.not. CellIsEmpty)then
                        call createIons(Nx,j,n,q,EdgeListHead)
                    end if
                end do
            end do

            q => EdgeListHead
            p => EdgeListHead%next
            do while(associated(p))
                p%x = p%x+p%vx*dt
                p%y = p%y+p%vy*dt
                if(p%x<xMaxBox .and. p%y>yMinBox .and. p%y<yMaxBox)then
                    q => p
                    p => p%next
                else
                    call deleteIon(EdgeListHead,q)
                    p => q%next
                end if
            end do
            call moveListA2B(EdgeListHead,IonListHead)
        end if

        if(yMinBoundary==4 .and. ppj==0)then
            call deleteList(MPIListHeadDown)
            q => EdgeListHead
            do n=1,NSpecies
                do i=0,Nx-1
                    CellIsEmpty = isCellEmpty(i,-1,n)
                    if(.not. CellIsEmpty)then
                        call createIons(i,-1,n,q,EdgeListHead)
                    end if
                end do
            end do

            q => EdgeListHead
            p => EdgeListHead%next
            do while(associated(p))
                p%x = p%x+p%vx*dt
                p%y = p%y+p%vy*dt
                if(p%y>yMinBox .and. p%x>xMinBox .and. p%x<xMaxBox)then
                    q => p
                    p => p%next
                else
                    call deleteIon(EdgeListHead,q)
                    p => q%next
                end if
            end do
            call moveListA2B(EdgeListHead,IonListHead)
        end if

        if(yMaxBoundary==4 .and. ppj==ppy-1)then
            call deleteList(MPIListHeadUp)
            q => EdgeListHead
            do n=1,NSpecies
                do i=0,Nx-1
                    CellIsEmpty = isCellEmpty(i,Ny,n)
                    if(.not. CellIsEmpty)then
                        call createIons(i,Ny,n,q,EdgeListHead)
                    end if
                end do
            end do

            q => EdgeListHead
            p => EdgeListHead%next
            do while(associated(p))
                p%x = p%x+p%vx*dt
                p%y = p%y+p%vy*dt
                if(p%y<yMaxBox .and. p%x>xMinBox .and. p%x<xMaxBox)then
                    q => p
                    p => p%next
                else
                    call deleteIon(EdgeListHead,q)
                    p => q%next
                end if
            end do
            call moveListA2B(EdgeListHead,IonListHead)
        end if

        ! Corner Ions
        if(xMinBoundary==4 .or. yMinBoundary==4)then
            if(ppi==0 .and. ppj==0)then
                q => CornerListHead
                do n=1,NSpecies
                    CellIsEmpty = isCellEmpty(-1,-1,n)
                    if(.not. CellIsEmpty)then
                        call createIons(-1,-1,n,q,CornerListHead)
                    end if
                end do
    
                q => CornerListHead
                p => CornerListHead%next
                do while(associated(p))
                    p%x = p%x+p%vx*dt
                    p%y = p%y+p%vy*dt
                    if(p%x>xMinBox .and. p%y>yMinBox)then
                        q => p
                        p => p%next
                    else
                        call deleteIon(CornerListHead,q)
                        p => q%next
                    end if
                end do
                call moveListA2B(CornerListHead,IonListHead)
            end if
        end if

        if(xMinBoundary==4 .or. yMaxBoundary==4)then
            if(ppi==0 .and. ppj==ppy-1)then
                q => CornerListHead
                do n=1,NSpecies
                    CellIsEmpty = isCellEmpty(-1,Ny,n)
                    if(.not. CellIsEmpty)then
                        call createIons(-1,Ny,n,q,CornerListHead)
                    end if
                end do
    
                q => CornerListHead
                p => CornerListHead%next
                do while(associated(p))
                    p%x = p%x+p%vx*dt
                    p%y = p%y+p%vy*dt
                    if(p%x>xMinBox .and. p%y<yMaxBox)then
                        q => p
                        p => p%next
                    else
                        call deleteIon(CornerListHead,q)
                        p => q%next
                    end if
                end do
                call moveListA2B(CornerListHead,IonListHead)
            end if
        end if

        if(xMaxBoundary==4 .or. yMinBoundary==4)then
            if(ppi==ppx-1 .and. ppj==0)then
                q => CornerListHead
                do n=1,NSpecies
                    CellIsEmpty = isCellEmpty(Nx,-1,n)
                    if(.not. CellIsEmpty)then
                        call createIons(Nx,-1,n,q,CornerListHead)
                    end if
                end do
    
                q => CornerListHead
                p => CornerListHead%next
                do while(associated(p))
                    p%x = p%x+p%vx*dt
                    p%y = p%y+p%vy*dt
                    if(p%x<xMaxBox .and. p%y>yMinBox)then
                        q => p
                        p => p%next
                    else
                        call deleteIon(CornerListHead,q)
                        p => q%next
                    end if
                end do
                call moveListA2B(CornerListHead,IonListHead)
            end if
        end if

        if(xMaxBoundary==4 .or. yMaxBoundary==4)then
            if(ppi==ppx-1 .and. ppj==ppy-1)then
                q => CornerListHead
                do n=1,NSpecies
                    CellIsEmpty = isCellEmpty(Nx,Ny,n)
                    if(.not. CellIsEmpty)then
                        call createIons(Nx,Ny,n,q,CornerListHead)
                    end if
                end do
    
                q => CornerListHead
                p => CornerListHead%next
                do while(associated(p))
                    p%x = p%x+p%vx*dt
                    p%y = p%y+p%vy*dt
                    if(p%x<xMaxBox .and. p%y<yMaxBox)then
                        q => p
                        p => p%next
                    else
                        call deleteIon(CornerListHead,q)
                        p => q%next
                    end if
                end do
                call moveListA2B(CornerListHead,IonListHead)
            end if
        end if

        deallocate(EdgeListHead)
        deallocate(CornerListHead)
    end subroutine applyThermalForParticle

    ! Useful Functions
    subroutine setPeriodicPosition(ListHead)
        implicit none
        type(PseudoParticle), pointer :: ListHead
        type(PseudoParticle), pointer :: p, q

        if(associated(ListHead,MPIListHeadLeft) .and. ppi==0)then
            q => ListHead
            p => ListHead%next
            do while(associated(p))
                p%x = p%x+(xMaxBox-xMinBox)
                q => p
                p => p%next
            end do
        end if

        if(associated(ListHead,MPIListHeadRight) .and. ppi==ppx-1)then
            q => ListHead
            p => ListHead%next
            do while(associated(p))
                p%x = p%x-(xMaxBox-xMinBox)
                q => p
                p => p%next
            end do
        end if

        if(associated(ListHead,MPIListHeadUp) .and. ppj==ppy-1)then
            q => ListHead
            p => ListHead%next
            do while(associated(p))
                p%y = p%y-(yMaxBox-yMinBox)
                q => p
                p => p%next
            end do
        end if

        if(associated(ListHead,MPIListHeadDown) .and. ppj==0)then
            q => ListHead
            p => ListHead%next
            do while(associated(p))
                p%y = p%y+(yMaxBox-yMinBox)
                q => p
                p => p%next
            end do
        end if
    end subroutine setPeriodicPosition

    subroutine useOpen3ForGrid(opt)
        implicit none
        integer(kind=4) :: opt, i

        if(xMinBoundary==opt)then
            if(ppi==0)then
                ux(0,:) = ux(0,:)*2.0d0
                uy(0,:) = uy(0,:)*2.0d0
                uz(0,:) = uz(0,:)*2.0d0
                ! ni(0,:) = ni(0,:)*2.0d0
                ! ne(0,:) = ne(0,:)*2.0d0
                ni(0,:) = ni(0,:)+ne0(0,:)/Z_avg/2
                ne(0,:) = ne(0,:)+ne0(0,:)/2
                if(ppj==ppy-1)then
                    ni(0,Ny) = ni(0,Ny)-ne0(0,Ny)/4/Z_avg
                    ne(0,Ny) = ne(0,Ny)-ne0(0,Ny)/4
                elseif(ppj==0)then
                    ni(0,0) = ni(0,0)-ne0(0,0)/4/Z_avg
                    ne(0,0) = ne(0,0)-ne0(0,0)/4
                end if
                do i=-1-ncx,-1
                    ux(i,:) = ux(0,:)
                    uy(i,:) = uy(0,:)
                    uz(i,:) = uz(0,:)
                    ni(i,:) = ne0(0,:)/Z_avg
                    ne(i,:) = ne0(0,:)
                    ! Bx(i,:) = 0.0d0
                    ! By(i,:) = 0.0d0
                    ! Bz(i,:) = 0.0d0
                    Bx(i,:) = Bx(0,:)
                    By(i,:) = By(1,:)
                    Bz(i,:) = Bz(1,:)
                    ! Ex(i,:) = Ex(1,:)
                    ! Ey(i,:) = Ey(0,:)
                    ! Ez(i,:) = Ez(0,:)
                    ! Ex(i,:) = 0.0d0
                    ! Ey(i,:) = 0.0d0
                    ! Ez(i,:) = 0.0d0
                end do
                ! Bx(-1-ncx:-1,:) = Bx0(-1-ncx:-1,:)
                ! By(-1-ncx:-1,:) = By0(-1-ncx:-1,:)
                ! Bz(-1-ncx:-1,:) = Bz0(-1-ncx:-1,:)
            end if
        end if

        if(xMaxBoundary==opt)then
            if(ppi==ppx-1)then
                ux(Nx,:) = ux(Nx,:)*2.0d0
                uy(Nx,:) = uy(Nx,:)*2.0d0
                uz(Nx,:) = uz(Nx,:)*2.0d0
                ! ni(Nx,:) = ni(Nx,:)*2.0d0
                ! ne(Nx,:) = ne(Nx,:)*2.0d0
                ni(Nx,:) = ni(Nx,:)+ne0(Nx,:)/Z_avg/2
                ne(Nx,:) = ne(Nx,:)+ne0(Nx,:)/2
                if(ppj==ppy-1)then
                    ni(Nx,Ny) = ni(Nx,Ny)-ne0(Nx,Ny)/4/Z_avg
                    ne(Nx,Ny) = ne(Nx,Ny)-ne0(Nx,Ny)/4
                elseif(ppj==0)then
                    ni(Nx,0) = ni(Nx,0)-ne0(Nx,0)/4/Z_avg
                    ne(Nx,0) = ne(Nx,0)-ne0(Nx,0)/4
                end if
                do i=Nx+1,Nx+ncx+1
                    ux(i,:) = ux(Nx,:)
                    uy(i,:) = uy(Nx,:)
                    uz(i,:) = uz(Nx,:)
                    ni(i,:) = ne0(Nx,:)/Z_avg
                    ne(i,:) = ne0(Nx,:)
                    ! Bx(i,:) = 0.0d0
                    ! By(i,:) = 0.0d0
                    ! Bz(i,:) = 0.0d0
                    Bx(i,:) = Bx(Nx,:)
                    By(i,:) = By(Nx,:)
                    Bz(i,:) = Bz(Nx,:)
                end do
            end if
        end if

        if(yMinBoundary==opt)then
            if(ppj==0)then
                ux(:,0) = ux(:,0)*2.0d0
                uy(:,0) = uy(:,0)*2.0d0
                uz(:,0) = uz(:,0)*2.0d0
                ! ni(:,0) = ni(:,0)*2.0d0
                ! ne(:,0) = ne(:,0)*2.0d0
                ni(:,0) = ni(:,0)+ne0(:,0)/Z_avg/2
                ne(:,0) = ne(:,0)+ne0(:,0)/2
                if(ppi==ppx-1)then
                    ni(Nx,0) = ni(Nx,0)-ne0(Nx,0)/4/Z_avg
                    ne(Nx,0) = ne(Nx,0)-ne0(Nx,0)/4
                elseif(ppi==0)then
                    ni(0,0) = ni(0,0)-ne0(0,0)/4/Z_avg
                    ne(0,0) = ne(0,0)-ne0(0,0)/4
                end if
                do i=-1-ncy,-1
                    ux(:,i) = ux(:,0)
                    uy(:,i) = uy(:,0)
                    uz(:,i) = uz(:,0)
                    ni(:,i) = ne0(:,0)/Z_avg
                    ne(:,i) = ne0(:,0)
                    ! Bx(:,i) = 0.0d0
                    ! By(:,i) = 0.0d0
                    ! Bz(:,i) = 0.0d0
                    Bx(:,i) = Bx(:,1)
                    By(:,i) = By(:,0)
                    Bz(:,i) = Bz(:,1)
                end do
            end if
        end if

        if(yMaxBoundary==opt)then
            if(ppj==ppy-1)then
                ux(:,Ny) = ux(:,Ny)*2.0d0
                uy(:,Ny) = uy(:,Ny)*2.0d0
                uz(:,Ny) = uz(:,Ny)*2.0d0
                ! ni(:,Ny) = ni(:,Ny)*2.0d0
                ! ne(:,Ny) = ne(:,Ny)*2.0d0
                ni(:,Ny) = ni(:,Ny)+ne0(:,Ny)/Z_avg/2
                ne(:,Ny) = ne(:,Ny)+ne0(:,Ny)/2
                if(ppi==ppx-1)then
                    ni(Nx,Ny) = ni(Nx,Ny)-ne0(Nx,Ny)/4/Z_avg
                    ne(Nx,Ny) = ne(Nx,Ny)-ne0(Nx,Ny)/4
                elseif(ppi==0)then
                    ni(0,Ny) = ni(0,Ny)-ne0(0,Ny)/4/Z_avg
                    ne(0,Ny) = ne(0,Ny)-ne0(0,Ny)/4
                end if
                do i=Ny+1,Ny+ncy+1
                    ux(:,i) = ux(:,Ny)
                    uy(:,i) = uy(:,Ny)
                    uz(:,i) = uz(:,Ny)
                    ni(:,i) = ne0(:,Ny)/Z_avg
                    ne(:,i) = ne0(:,Ny)
                    ! Bx(:,i) = 0.0d0
                    ! By(:,i) = 0.0d0
                    ! Bz(:,i) = 0.0d0
                    Bx(:,i) = Bx(:,Ny)
                    By(:,i) = By(:,Ny)
                    Bz(:,i) = Bz(:,Ny)
                end do
            end if
        end if
        call useOpen3ForGrid_in_exchange(opt,Ex)
        call useOpen3ForGrid_in_exchange(opt,Ey)
        call useOpen3ForGrid_in_exchange(opt,Ez)
    end subroutine useOpen3ForGrid

    subroutine useReflectForGrid(opt)
        implicit none
        integer(kind=4) :: opt, i

        if(xMinBoundary==opt)then
            if(ppi==0)then
                ux(0,:) = 0.0d0
                uy(0,:) = uy(0,:)*2.0d0
                uz(0,:) = uz(0,:)*2.0d0
                ni(0,:) = ni(0,:)*2.0d0
                ne(0,:) = ne(0,:)*2.0d0
                do i=-1-ncx,-1
                    ux(i,:) = -ux(-i,:)
                    uy(i,:) = uy(-i,:)
                    uz(i,:) = uz(-i,:)
                    ni(i,:) = ni(-i,:)
                    ne(i,:) = ne(-i,:)
                    ! Ex(i,:) = -Ex(-i,:)
                    ! Ey(i,:) = Ey(-i,:)
                    ! Ez(i,:) = Ez(-i,:)
                    Bx(i,:) = Bx(-i,:)
                    By(i,:) = 2*By0(0,:)-By(-i,:)
                    Bz(i,:) = 2*Bz0(0,:)-Bz(-i,:)
                end do
            end if
        end if

        if(xMaxBoundary==opt)then
            if(ppi==ppx-1)then
                ux(Nx,:) = 0.0d0
                Bx(Nx,:) = 0.0d0
                uy(Nx,:) = uy(Nx,:)*2.0d0
                uz(Nx,:) = uz(Nx,:)*2.0d0
                ni(Nx,:) = ni(Nx,:)*2.0d0
                ne(Nx,:) = ne(Nx,:)*2.0d0
                do i=1,ncx+1
                    ux(Nx+i,:) = -ux(Nx-i,:)
                    uy(Nx+i,:) = uy(Nx-i,:)
                    uz(Nx+i,:) = uz(Nx-i,:)
                    ni(Nx+i,:) = ni(Nx-i,:)
                    ne(Nx+i,:) = ne(Nx-i,:)
                    ! Ex(Nx+i,:) = -Ex(Nx-i+1,:)
                    ! Ey(Nx+i,:) = -Ey(Nx-i,:)
                    ! Ez(Nx+i,:) = Ez(Nx-i,:)
                    Bx(Nx+i,:) = Bx(Nx-i,:)
                    By(Nx+i,:) = 2*By0(0,:)-By(Nx-i+1,:)
                    ! Bz(Nx+i,:) = 2*Bz0(0,:)-Bz(Nx-i+1,:)
                    Bz(Nx+i,:) = Bz(Nx-i+1,:)
                end do
            end if
        end if

        if(yMinBoundary==opt)then
            if(ppj==0)then
                ux(:,0) = ux(:,0)*2.0d0
                uy(:,0) = 0.0d0
                uz(:,0) = uz(:,0)*2.0d0
                ni(:,0) = ni(:,0)*2.0d0
                ne(:,0) = ne(:,0)*2.0d0
                do i=-1-ncy,-1
                    ux(:,i) = ux(:,-i)
                    uy(:,i) = -uy(:,-i)
                    uz(:,i) = uz(:,-i)
                    ni(:,i) = ni(:,-i)
                    ne(:,i) = ne(:,-i)
                    ! Ex(:,i) = Ex(:,-i)
                    ! Ey(:,i) = -Ey(:,-i)
                    ! Ez(:,i) = Ez(:,-i)
                    Bx(:,i) = 2*Bx0(0,:)-Bx(:,-i)
                    By(:,i) = By(:,-i)
                    Bz(:,i) = 2*Bz0(0,:)-Bz(:,-i)
                end do
            end if
        end if

        if(yMaxBoundary==opt)then
            if(ppj==ppy-1)then
                ux(:,Ny) = ux(:,Ny)*2.0d0
                uy(:,Ny) = 0.0d0
                uz(:,Ny) = uz(:,Ny)*2.0d0
                ni(:,Ny) = ni(:,Ny)*2.0d0
                ne(:,Ny) = ne(:,Ny)*2.0d0
                do i=1,ncy+1
                    ux(:,Ny+i) = ux(:,Ny-i)
                    uy(:,Ny+i) = -uy(:,Ny-i)
                    uz(:,Ny+i) = uz(:,Ny-i)
                    ni(:,Ny+i) = ni(:,Ny-i)
                    ne(:,Ny+i) = ne(:,Ny-i)
                    ! Ex(:,Ny+i) = Ex(:,Ny-i)
                    ! Ey(:,Ny+i) = -Ey(:,Ny-i+1)
                    ! Ez(:,Ny+i) = Ez(:,Ny-i)
                    Bx(:,Ny+i) = 2*Bx0(0,:)-Bx(:,Ny-i+1)
                    By(:,Ny+i) = By(:,Ny-i)
                    Bz(:,Ny+i) = 2*Bz0(0,:)-Bz(:,Ny-i+1)
                end do
            end if
        end if
        call useReflectForGrid_in_exchange(opt,Ex)
        call useReflectForGrid_in_exchange(opt,Ey)
        call useReflectForGrid_in_exchange(opt,Ez)
    end subroutine useReflectForGrid

    subroutine useXinverseForGrid(opt)
        implicit none
        integer(kind=4) :: opt, i

        if(xMinBoundary==opt)then
            if(ppi==0)then
                ux(0,:) = 0.0d0
                uy(0,:) = uy(0,:)*2.0d0
                uz(0,:) = uz(0,:)*2.0d0
                ni(0,:) = ni(0,:)*2.0d0
                ne(0,:) = ne(0,:)*2.0d0
                ! Ex(0,:) = 0.0d0
                ! Ey(0,:) = 0.0d0
                
                By(0,:) = 0.0d0
                Bz(0,:) = 0.0d0
                do i=-1-ncx,-1
                    ux(i,:) = -ux(-i,:)
                    uy(i,:) = uy(-i,:)
                    uz(i,:) = uz(-i,:)
                    ni(i,:) = ni(-i,:)
                    ne(i,:) = ne(-i,:)
                    Ex(i,:) = -Ex(-i,:)
                    Ey(i,:) = Ey(-i,:)
                    Ez(i,:) = Ez(-i,:)
                    Bx(i,:) = Bx(-i,:)
                    By(i,:) = -By(-i,:)
                    Bz(i,:) = -Bz(-i,:)
                end do
            end if
        end if

        if(xMaxBoundary==opt)then
            if(ppi==ppx-1)then
                ux(Nx,:) = 0.0d0
                uy(Nx,:) = uy(Nx,:)*2.0d0
                uz(Nx,:) = uz(Nx,:)*2.0d0
                ni(Nx,:) = ni(Nx,:)*2.0d0
                ne(Nx,:) = ne(Nx,:)*2.0d0
                ! Ex(Nx,:) = 0.0d0
                ! Ey(Nx,:) = 0.0d0
                
                By(Nx,:) = 0.0d0
                Bz(Nx,:) = 0.0d0
                do i=1,ncx+1
                    ux(Nx+i,:) = -ux(Nx-i,:)
                    uy(Nx+i,:) = uy(Nx-i,:)
                    uz(Nx+i,:) = uz(Nx-i,:)
                    ni(Nx+i,:) = ni(Nx-i,:)
                    ne(Nx+i,:) = ne(Nx-i,:)
                    Ex(Nx+i,:) = -Ex(Nx-i+1,:)
                    Ey(Nx+i,:) = Ey(Nx-i,:)
                    Ez(Nx+i,:) = Ez(Nx-i,:)
                    Bx(Nx+i,:) = Bx(Nx-i,:)
                    By(Nx+i,:) = -By(Nx-i+1,:)
                    Bz(Nx+i,:) = -Bz(Nx-i+1,:)
                end do
            end if
        end if

        if(yMinBoundary==opt)then
            if(ppj==0)then
                ux(:,0) = ux(:,0)*2.0d0
                uy(:,0) = 0.0d0
                uz(:,0) = uz(:,0)*2.0d0
                ni(:,0) = ni(:,0)*2.0d0
                ne(:,0) = ne(:,0)*2.0d0
                ! Ex(:,0) = 0.0d0
                ! Ey(:,0) = 0.0d0
                
                Bx(:,0) = 0.0d0
                Bz(:,0) = 0.0d0
                do i=-1-ncy,-1
                    ux(:,i) = ux(:,-i)
                    uy(:,i) = -uy(:,-i)
                    uz(:,i) = uz(:,-i)
                    ni(:,i) = ni(:,-i)
                    ne(:,i) = ne(:,-i)
                    Ex(:,i) = Ex(:,-i)
                    Ey(:,i) = -Ey(:,-i)
                    Ez(:,i) = Ez(:,-i)
                    Bx(:,i) = -Bx(:,-i)
                    By(:,i) = By(:,-i)
                    Bz(:,i) = -Bz(:,-i)
                end do
            end if
        end if

        if(yMaxBoundary==opt)then
            if(ppj==ppy-1)then
                ux(:,Ny) = ux(:,Ny)*2.0d0
                uy(:,Ny) = 0.0d0
                uz(:,Ny) = uz(:,Ny)*2.0d0
                ni(:,Ny) = ni(:,Ny)*2.0d0
                ne(:,Ny) = ne(:,Ny)*2.0d0
                ! Ex(:,Ny) = 0.0d0
                ! Ey(:,Ny) = 0.0d0
                
                Bx(:,Ny) = 0.0d0
                Bz(:,Ny) = 0.0d0
                do i=1,ncy+1
                    ux(:,Ny+i) = ux(:,Ny-i)
                    uy(:,Ny+i) = -uy(:,Ny-i)
                    uz(:,Ny+i) = uz(:,Ny-i)
                    ni(:,Ny+i) = ni(:,Ny-i)
                    ne(:,Ny+i) = ne(:,Ny-i)
                    Ex(:,Ny+i) = Ex(:,Ny-i)
                    Ey(:,Ny+i) = -Ey(:,Ny-i+1)
                    Ez(:,Ny+i) = Ez(:,Ny-i)
                    Bx(:,Ny+i) = -Bx(:,Ny-i+1)
                    By(:,Ny+i) = By(:,Ny-i)
                    Bz(:,Ny+i) = -Bz(:,Ny-i+1)
                end do
            end if
        end if
    end subroutine useXinverseForGrid

    subroutine useOpen3ForGrid_in_exchange(opt,G)
        implicit none
        integer(kind=4) :: opt, i
        real(kind=8) :: G(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)

        if(xMinBoundary==opt)then
            if(ppi==0)then
                if(loc(G)==loc(By) .or. loc(G)==loc(Bz) .or. &
                loc(G)==loc(Bfy) .or. loc(G)==loc(Bfz) .or. &
                loc(G)==loc(Bwy) .or. loc(G)==loc(Bwz))then
                    do i=-1-ncx,-1
                        G(i,:) = G(1,:)
                        
                    end do
                    ! if(loc(G)==loc(Bx) .or.loc(G)==loc(Bfx)) G(-1-ncx:-1,:) = Bx0(-1-ncx:-1,:)
                    ! if(loc(G)==loc(By) .or.loc(G)==loc(Bfy)) G(-1-ncx:-1,:) = By0(-1-ncx:-1,:)
                    ! if(loc(G)==loc(Bz) .or.loc(G)==loc(Bfz)) G(-1-ncx:-1,:) = Bz0(-1-ncx:-1,:)

                elseif(loc(G)==loc(Ex) .or. loc(G)==loc(Ey) .or. loc(G)==loc(Ez))then
                    ! G(-1-ncx:-1,:) = 0.0d0
                    call solveOhmEquationfor_Xboundary(0,G) ! 0 means MinBoundary, 1 means MaxBoundary.
                    G(-1-ncx:-1,Ny+1:Ny+ncy+1) = 0.0d0
                    G(-1-ncx:-1,-1-ncy:-1) = 0.0d0
                    ! call applyGridCommunication(G)
                elseif(loc(G)==loc(ni))then
                    do i=-1-ncx,-1
                        G(i,:) = ne0(0,:)/Z_avg
                        
                    end do
                elseif(loc(G)==loc(ne))then
                    do i=-1-ncx,-1
                        G(i,:) = ne0(0,:)
                        
                    end do
                
                else
                    do i=-1-ncx,-1
                        G(i,:) = G(0,:)
                        
                    end do
                end if
                
            end if
        end if

        if(xMaxBoundary==opt)then
            if(ppi==ppx-1)then
                
                if(loc(G)==loc(Bx) .or. loc(G)==loc(By) .or. loc(G)==loc(Bz) .or. &
                loc(G)==loc(Bfx) .or. loc(G)==loc(Bfy) .or. loc(G)==loc(Bfz))then
                    
                    if(loc(G)==loc(Bx) .or.loc(G)==loc(Bfx)) G(Nx+1:Nx+ncx+1,:) = Bx0(Nx+1:Nx+ncx+1,:)
                    if(loc(G)==loc(By) .or.loc(G)==loc(Bfy)) G(Nx+1:Nx+ncx+1,:) = By0(Nx+1:Nx+ncx+1,:)
                    if(loc(G)==loc(Bz) .or.loc(G)==loc(Bfz)) G(Nx+1:Nx+ncx+1,:) = Bz0(Nx+1:Nx+ncx+1,:)
                    
                elseif(loc(G)==loc(Ex) .or. loc(G)==loc(Ey) .or. loc(G)==loc(Ez))then
                    call solveOhmEquationfor_Xboundary(1,G)
                    G(Nx+1:Nx+ncx+1,Ny+1:Ny+ncy+1) = 0.0d0
                    G(Nx+1:Nx+ncx+1,-1-ncy:-1) = 0.0d0
                    ! call applyGridCommunication(G)
                elseif(loc(G)==loc(ni))then
                    do i=Nx+1,Nx+ncx+1
                        G(i,:) = ne0(Nx,:)/Z_avg
                        
                    end do
                elseif(loc(G)==loc(ne))then
                    do i=Nx+1,Nx+ncx+1
                        G(i,:) = ne0(Nx,:)
                        
                    end do
                else
                    do i=Nx+1,Nx+ncx+1
                        G(i,:) = G(Nx,:)
                    end do
                end if
            end if
        end if

        if(yMinBoundary==opt)then
            if(ppj==0)then

                if(loc(G)==loc(Bx) .or. loc(G)==loc(By) .or. loc(G)==loc(Bz) .or. &
                loc(G)==loc(Bfx) .or. loc(G)==loc(Bfy) .or. loc(G)==loc(Bfz))then
                    
                    if(loc(G)==loc(Bx) .or.loc(G)==loc(Bfx)) G(:,-1-ncy:-1) = Bx0(:,-1-ncy:-1)
                    if(loc(G)==loc(By) .or.loc(G)==loc(Bfy)) G(:,-1-ncy:-1) = By0(:,-1-ncy:-1)
                    if(loc(G)==loc(Bz) .or.loc(G)==loc(Bfz)) G(:,-1-ncy:-1) = Bz0(:,-1-ncy:-1)
                    
                elseif(loc(G)==loc(Ex) .or. loc(G)==loc(Ey) .or. loc(G)==loc(Ez))then
                    call solveOhmEquationfor_Yboundary(0,G)
                    G(Nx+1:Nx+ncx+1,-1-ncy:-1) = 0.0d0
                    G(-1-ncx:-1,-1-ncy:-1) = 0.0d0
                    ! call applyGridCommunication(G)
                elseif(loc(G)==loc(ni))then
                    do i=-1-ncy,-1
                        G(:,i) = ne0(:,0)/Z_avg
                        
                    end do
                elseif(loc(G)==loc(ne))then
                    do i=-1-ncy,-1
                        G(:,i) = ne0(:,0)
                        
                    end do
                else
                    do i=-1-ncy,-1
                        G(:,i) = G(:,0)
                    end do
                end if

            end if
        end if

        if(yMaxBoundary==opt)then
            if(ppj==ppy-1)then
                
                if(loc(G)==loc(Bx) .or. loc(G)==loc(By) .or. loc(G)==loc(Bz) .or. &
                loc(G)==loc(Bfx) .or. loc(G)==loc(Bfy) .or. loc(G)==loc(Bfz))then
                    
                    if(loc(G)==loc(Bx) .or.loc(G)==loc(Bfx)) G(:,Ny+1:Ny+ncy+1) = Bx0(:,Ny+1:Ny+ncy+1)
                    if(loc(G)==loc(By) .or.loc(G)==loc(Bfy)) G(:,Ny+1:Ny+ncy+1) = By0(:,Ny+1:Ny+ncy+1)
                    if(loc(G)==loc(Bz) .or.loc(G)==loc(Bfz)) G(:,Ny+1:Ny+ncy+1) = Bz0(:,Ny+1:Ny+ncy+1)
                    
                elseif(loc(G)==loc(Ex) .or. loc(G)==loc(Ey) .or. loc(G)==loc(Ez))then
                    
                    call solveOhmEquationfor_Yboundary(1,G)
                    
                    G(Nx+1:Nx+ncx+1,Ny+1:Ny+ncy+1) = 0.0d0
                    G(-1-ncx:-1,Ny+1:Ny+ncy+1) = 0.0d0
                    ! call applyGridCommunication(G)
                elseif(loc(G)==loc(ni))then
                    do i=Ny+1,Ny+ncy+1
                        G(:,i) = ne0(:,Ny)/Z_avg
                        
                    end do
                elseif(loc(G)==loc(ne))then
                    do i=Ny+1,Ny+ncy+1
                        G(:,i) = ne0(:,Ny)
                        
                    end do
                else
                    do i=Ny+1,Ny+ncy+1
                        G(:,i) = G(:,Ny)
                    end do
                end if

            end if
        end if
    end subroutine useOpen3ForGrid_in_exchange

    subroutine solveOhmEquationfor_Xboundary(opt,G)
        implicit none
        integer(kind=4) :: opt, i
        real(kind=8) :: G(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)
        real(kind=8) :: p_ele(1:ncx+1,-1-ncy:Ny+ncy+1)
        real(kind=8) :: Bfield1(1:ncx+1,-ncy:Ny+ncy+1)
        real(kind=8) :: Bfield2(1:ncx+1,-1-ncy:Ny+ncy+1)
        real(kind=8) :: Bfield3(1:ncx+1,-ncy:Ny+ncy+1)
        real(kind=8) :: temp_bx(1+ncx,-ncy:Ny+ncy)
        real(kind=8) :: temp_by(1+ncx,-ncy:Ny+ncy)
        real(kind=8) :: temp_bz(1+ncx,-ncy:Ny+ncy)
        real(kind=8) :: bgx(-1-ncx:Nx+ncx+1,-ncy:Ny+ncy+1), bgy(-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)
        real(kind=8) :: bgz(-ncx:Nx+ncx+1,-ncy:Ny+ncy+1)
        real(kind=8) :: velx(1:ncx+1,-1-ncy:Ny+ncy+1), uex(1:ncx+1,-1-ncy:Ny+ncy+1)
        real(kind=8) :: vely(1:ncx+1,-1-ncy:Ny+ncy+1), uey(1:ncx+1,-1-ncy:Ny+ncy+1)
        real(kind=8) :: velz(1:ncx+1,-1-ncy:Ny+ncy+1), uez(1:ncx+1,-1-ncy:Ny+ncy+1)
        real(kind=8) :: n_ele(1:ncx+1,-1-ncy:Ny+ncy+1)
        real(kind=8) :: n_ion(1:ncx+1,-1-ncy:Ny+ncy+1)
        real(kind=8) :: Temp(ncx+1,-ncy:Ny+ncy,3)
        
        real(kind=8) :: Temp1(ncx+1,Ny+1), Temp2(ncx+1,Ny+1), Temp3(ncx+1,Ny+1)
        real(kind=8) :: curl(ncx+1,Ny+1,3)
        
        bgx = getBxGrid(Bx)
        bgy = getByGrid(By)
        bgz = getBzGrid(Bz)
        if(opt==0)then
            p_ele = pe(-1-ncx:-1,:)
            n_ele = ne(-1-ncx:-1,:)
            n_ion = ni(-1-ncx:-1,:)
            velx = ux(-1-ncx:-1,:)
            vely = uy(-1-ncx:-1,:)
            velz = uz(-1-ncx:-1,:)
            Bfield1 = bgx(-1-ncx:-1,:)
            Bfield2 = bgy(-ncx:0,:)
            Bfield3 = bgz(-ncx:0,:)
        else
            p_ele = pe(Nx+1:Nx+ncx+1,:)
            n_ele = ne(Nx+1:Nx+ncx+1,:)
            n_ion = ni(Nx+1:Nx+ncx+1,:)
            velx = ux(Nx+1:Nx+ncx+1,:)
            vely = uy(Nx+1:Nx+ncx+1,:)
            velz = uz(Nx+1:Nx+ncx+1,:)
            Bfield1 = bgx(Nx+1:Nx+ncx+1,:)
            Bfield2 = bgy(Nx+1:Nx+ncx+1,:)
            Bfield3 = bgz(Nx+1:Nx+ncx+1,:)
        end if
        temp_bx = (Bfield1(:,-ncy:Ny+ncy)+Bfield1(:,1-ncy:Ny+ncy+1))/2.0
        temp_by = Bfield2(:,-ncy:Ny+ncy)
        temp_bz = (Bfield3(:,-ncy:Ny+ncy)+Bfield3(:,1-ncy:Ny+ncy+1))/2.0
        
        Temp = cross_product_arr2(velx(:,-ncy:Ny+ncy),vely(:,-ncy:Ny+ncy),velz(:,-ncy:Ny+ncy) &
            ,temp_bx,temp_by,temp_bz,ncx+1,Ny+2*ncy+1)
    
        Temp1 = (Bfield3(:,1:Ny+1)-Bfield3(:,0:Ny))/dy! DyBz-DzBy
        Temp2 = 0.0d0
        Temp3 = -(Bfield1(:,1:Ny+1)-Bfield1(:,0:Ny))/dy! DxBy-DyBx
        
        curl = cross_product_arr2(Temp1,Temp2,Temp3, &
                temp_bx(:,0:Ny),temp_by(:,0:Ny),temp_bz(:,0:Ny),ncx+1,Ny+1)
    
        if(opt==0)then
            if(loc(G)==loc(Ex))then
                G(-1-ncx:-1,0:Ny) = - Temp(:,0:Ny,1) + curl(:,:,1)/n_ion(:,0:Ny)/Z_avg
            elseif(loc(G)==loc(Ey))then
                G(-1-ncx:-1,1:Ny) = 2*(p_ele(:,0:Ny-1)-p_ele(:,1:Ny))/(n_ele(:,0:Ny-1)+n_ele(:,1:Ny))/dy &
                            - (Temp(:,1:Ny,2)+Temp(:,0:Ny-1,2))/2 &
                            + (curl(:,1:Ny,2)+curl(:,2:Ny+1,2))/(n_ion(:,0:Ny-1)+n_ion(:,1:Ny))/Z_avg
            elseif(loc(G)==loc(Ez))then
                G(-1-ncx:-1,0:Ny) = -Temp(:,0:Ny,3) + curl(:,:,3)/n_ion(:,0:Ny)/Z_avg
            end if
        else
            if(loc(G)==loc(Ex))then
                G(Nx+1:Nx+ncx+1,0:Ny) = - Temp(:,0:Ny,1) + curl(:,:,1)/n_ion(:,0:Ny)/Z_avg
            elseif(loc(G)==loc(Ey))then
                G(Nx+1:Nx+ncx+1,1:Ny) = 2*(p_ele(:,0:Ny-1)-p_ele(:,1:Ny))/(n_ele(:,0:Ny-1)+n_ele(:,1:Ny))/dy &
                - (Temp(:,1:Ny,2)+Temp(:,0:Ny-1,2))/2 &
                + (curl(:,1:Ny,2)+curl(:,2:Ny+1,2))/(n_ion(:,0:Ny-1)+n_ion(:,1:Ny))/Z_avg
            elseif(loc(G)==loc(Ez))then
                G(Nx+1:Nx+ncx+1,0:Ny) = -Temp(:,0:Ny,3) + curl(:,:,3)/n_ion(:,0:Ny)/Z_avg
            end if
        end if

    end subroutine solveOhmEquationfor_Xboundary

    subroutine solveOhmEquationfor_Yboundary(opt,G)
        implicit none
        integer(kind=4) :: opt, i
        real(kind=8) :: G(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)
        real(kind=8) :: p_ele(-1-ncx:Nx+ncx+1,ncy+1)
        real(kind=8) :: Bfield1(-1-ncx:Nx+ncx+1,ncy+1)
        real(kind=8) :: Bfield2(-ncx:Nx+ncx+1,ncy+1)
        real(kind=8) :: Bfield3(-ncx:Nx+ncx+1,ncy+1)
        real(kind=8) :: temp_bx(-ncx:Nx+ncx,ncy+1)
        real(kind=8) :: temp_by(-ncx:Nx+ncx,ncy+1)
        real(kind=8) :: temp_bz(-ncx:Nx+ncx,ncy+1)
        real(kind=8) :: bgx(-1-ncx:Nx+ncx+1,-ncy:Ny+ncy+1), bgy(-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)
        real(kind=8) :: bgz(-ncx:Nx+ncx+1,-ncy:Ny+ncy+1)
        real(kind=8) :: velx(-1-ncx:Nx+ncx+1,ncy+1)
        real(kind=8) :: vely(-1-ncx:Nx+ncx+1,ncy+1)
        real(kind=8) :: velz(-1-ncx:Nx+ncx+1,ncy+1)
        real(kind=8) :: n_ele(-1-ncx:Nx+ncx+1,ncy+1)
        real(kind=8) :: n_ion(-1-ncx:Nx+ncx+1,ncy+1)
        real(kind=8) :: Temp(-ncx:Nx+ncx,ncy+1,3)
        
        real(kind=8) :: Temp1(Nx+1,ncy+1), Temp2(Nx+1,ncy+1), Temp3(Nx+1,ncy+1)
        real(kind=8) :: curl(Nx+1,ncy+1,3)
        
        bgx = getBxGrid(Bx)
        bgy = getByGrid(By)
        bgz = getBzGrid(Bz)
        if(opt==0)then
            p_ele = pe(:,-1-ncy:-1)
            n_ele = ne(:,-1-ncy:-1)
            n_ion = ni(:,-1-ncy:-1)
            velx = ux(:,-1-ncy:-1)
            vely = uy(:,-1-ncy:-1)
            velz = uz(:,-1-ncy:-1)
            Bfield1 = bgx(:,-ncy:0)
            Bfield2 = bgy(:,-1-ncy:-1)
            Bfield3 = bgz(:,-ncy:0)
        else
            p_ele = pe(:,Ny+1:Ny+ncy+1)
            n_ele = ne(:,Ny+1:Ny+ncy+1)
            n_ion = ni(:,Ny+1:Ny+ncy+1)
            velx = ux(:,Ny+1:Ny+ncy+1)
            vely = uy(:,Ny+1:Ny+ncy+1)
            velz = uz(:,Ny+1:Ny+ncy+1)
            Bfield1 = bgx(:,Ny+1:Ny+ncy+1)
            Bfield2 = bgy(:,Ny+1:Ny+ncy+1)
            Bfield3 = bgz(:,Ny+1:Ny+ncy+1)
        end if
    
        temp_bx = Bfield1(-ncx:Nx+ncx,:)
        temp_by = (Bfield2(-ncx:Nx+ncx,:)+Bfield2(1-ncx:Nx+ncx+1,:))/2.0
        temp_bz = (Bfield3(-ncx:Nx+ncx,:)+Bfield3(1-ncx:Nx+ncx+1,:))/2.0
    
        Temp = cross_product_arr2(velx(-ncx:nx+ncx,:),vely(-ncx:nx+ncx,:),velz(-ncx:nx+ncx,:) &
            ,temp_bx,temp_by,temp_bz,Nx+2*ncx+1,ncy+1)
    
        Temp1 = 0.0d0! DyBz-DzBy
        Temp2 = (Bfield3(0:Nx,:)-Bfield3(1:Nx+1,:))/2/dx! DzBx-DxBz
        Temp3 = (Bfield2(1:Nx+1,:)-Bfield2(0:Nx,:))/dx! DxBy-DyBx
        
        curl = cross_product_arr2(Temp1,Temp2,Temp3, &
                temp_bx(0:Nx,:),temp_by(0:Nx,:),temp_bz(0:Nx,:),Nx+1,ncy+1)
    
        if(opt==0)then
            if(loc(G)==loc(Ex))then
                G(1:Nx,-1-ncy:-1) = 2*(p_ele(0:Nx-1,:)-p_ele(1:Nx,:))/(n_ele(0:Nx-1,:)+n_ele(1:Nx,:))/dx &
                - (Temp(0:Nx-1,:,1)+Temp(1:Nx,:,1))/2 + (curl(1:Nx,:,1)+curl(2:Nx+1,:,1))/(n_ion(0:Nx-1,:)+n_ion(1:Nx,:))/Z_avg
            elseif(loc(G)==loc(Ey))then
                G(0:Nx,-1-ncy:-1) = - Temp(0:Nx,:,2) + curl(:,:,2)/n_ion(0:Nx,:)/Z_avg
            elseif(loc(G)==loc(Ez))then
                G(0:Nx,-1-ncy:-1) = -Temp(0:Nx,:,3) + curl(:,:,3)/n_ion(0:Nx,:)/Z_avg
            end if
        else
            if(loc(G)==loc(Ex))then
                G(0:Nx,Ny+1:Ny+ncy+1) = 2*(p_ele(0:Nx-1,:)-p_ele(1:Nx,:))/(n_ele(0:Nx-1,:)+n_ele(1:Nx,:))/dx &
                - (Temp(0:Nx-1,:,1)+Temp(1:Nx,:,1))/2 + (curl(1:Nx,:,1)+curl(2:Nx+1,:,1))/(n_ion(0:Nx-1,:)+n_ion(1:Nx,:))/Z_avg
            elseif(loc(G)==loc(Ey))then
                G(0:Nx,Ny+1:Ny+ncy+1) = - Temp(0:Nx,:,2) + curl(:,:,2)/n_ion(0:Nx,:)/Z_avg
            elseif(loc(G)==loc(Ez))then
                G(0:Nx,Ny+1:Ny+ncy+1) = -Temp(0:Nx,:,3) + curl(:,:,3)/n_ion(0:Nx,:)/Z_avg
            end if
        end if
        
    end subroutine solveOhmEquationfor_Yboundary

    subroutine useReflectForGrid_in_exchange(opt,G)
        implicit none
        integer(kind=4) :: opt, i
        real(kind=8) :: G(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)
    
        if(xMinBoundary==opt)then
            if(ppi==0)then
                if(loc(G)==loc(By) .or. loc(G)==loc(Bz) .or. &
                loc(G)==loc(Bfy) .or. loc(G)==loc(Bfz) .or. &
                loc(G)==loc(Bwy) .or. loc(G)==loc(Bwz))then
                    do i=-1-ncx,-1
                        G(i,:) = -G(-i,:)
                    end do
                elseif(loc(G)==loc(Ex) .or. loc(G)==loc(Ey) .or. loc(G)==loc(Ez))then
                    call solveOhmEquationfor_Xboundary(0,G) ! 0 means MinBoundary, 1 means MaxBoundary.
                    G(-1-ncx:-1,Ny+1:Ny+ncy+1) = 0.0d0
                    G(-1-ncx:-1,-1-ncy:-1) = 0.0d0
                elseif(loc(G)==loc(ux))then
                    do i=-1-ncx,-1
                        G(i,:) = -G(-i,:)
                        
                    end do
                    G(0,:) = 0.0d0
                else
                    do i=-1-ncx,-1
                        G(i,:) = G(-i,:)
                        
                    end do
                end if
                
            end if
        end if
    
        if(xMaxBoundary==opt)then
            if(ppi==ppx-1)then
                
                if(loc(G)==loc(Bz) .or. &
                loc(G)==loc(Bfz) .or. &
                loc(G)==loc(Bwz))then
                    do i=Nx+1,Nx+ncx+1
                        G(i,:) = G(2*Nx-i+1,:)
                    end do
                elseif(loc(G)==loc(By) .or. &
                    loc(G)==loc(Bfy) .or. &
                    loc(G)==loc(Bwy))then
                        do i=Nx+1,Nx+ncx+1
                            G(i,:) = 2*By0(0,:)-G(2*Nx-i+1,:)
                        end do
                elseif(loc(G)==loc(Ex) .or. loc(G)==loc(Ey) .or. loc(G)==loc(Ez))then
                    call solveOhmEquationfor_Xboundary(1,G)
                    G(Nx+1:Nx+ncx+1,Ny+1:Ny+ncy+1) = 0.0d0
                    G(Nx+1:Nx+ncx+1,-1-ncy:-1) = 0.0d0
                elseif(loc(G)==loc(ux))then
                    do i=Nx+1,Nx+ncx+1
                        G(i,:) = -G(2*Nx-i,:)
                        
                    end do
                    G(Nx,:) = 0.0d0
                else
                    do i=Nx+1,Nx+ncx+1
                        G(i,:) = G(2*Nx-i,:)
                    end do
                end if
            end if
        end if
    
        if(yMinBoundary==opt)then
            if(ppj==0)then
    
                if(loc(G)==loc(Bx) .or. loc(G)==loc(Bz) .or. &
                loc(G)==loc(Bfx) .or. loc(G)==loc(Bfz) .or. &
                loc(G)==loc(Bwx) .or. loc(G)==loc(Bwz))then
                    do i=-1-ncy,-1
                        G(:,i) = -G(:,-i)
                    end do
                elseif(loc(G)==loc(Ex) .or. loc(G)==loc(Ey) .or. loc(G)==loc(Ez))then
                    call solveOhmEquationfor_Yboundary(0,G)
                    G(Nx+1:Nx+ncx+1,-1-ncy:-1) = 0.0d0
                    G(-1-ncx:-1,-1-ncy:-1) = 0.0d0
                elseif(loc(G)==loc(uy))then
                    do i=-1-ncy,-1
                        G(:,i) = -G(:,-i)
                    end do
                    G(:,0) = 0.0d0
                else
                    do i=-1-ncy,-1
                        G(:,i) = G(:,-i)
                        
                    end do
                end if
    
            end if
        end if
    
        if(yMaxBoundary==opt)then
            if(ppj==ppy-1)then
                
                
                if(loc(G)==loc(Bx) .or. loc(G)==loc(Bz) .or. &
                loc(G)==loc(Bfx) .or. loc(G)==loc(Bfz) .or. &
                loc(G)==loc(Bwx) .or. loc(G)==loc(Bwz))then
                    do i=Ny+1,Ny+ncy+1
                        G(:,i) = -G(:,2*Ny-i+1)
                    end do
                elseif(loc(G)==loc(Ex) .or. loc(G)==loc(Ey) .or. loc(G)==loc(Ez))then
                    
                    call solveOhmEquationfor_Yboundary(1,G)
                    
                    G(Nx+1:Nx+ncx+1,Ny+1:Ny+ncy+1) = 0.0d0
                    G(-1-ncx:-1,Ny+1:Ny+ncy+1) = 0.0d0
                elseif(loc(G)==loc(uy))then
                    do i=Ny+1,Ny+ncy+1
                        G(:,i) = -G(:,2*Ny-i)
                    end do
                    G(:,Ny) = 0.0d0
                else
                    do i=Ny+1,Ny+ncy+1
                        G(:,i) = G(:,2*Ny-i)
                    end do
                end if
    
            end if
        end if
    end subroutine useReflectForGrid_in_exchange

    subroutine useXinverseForGrid_in_exchange(opt,G)
        implicit none
        integer(kind=4) :: opt, i
        real(kind=8) :: G(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)

        if(xMinBoundary==opt)then
            if(ppi==0)then
                if(loc(G)==loc(Ex) .or. &
                   loc(G)==loc(By) .or. loc(G)==loc(Bz) .or. &
                   loc(G)==loc(Bfy) .or. loc(G)==loc(Bfz) .or. &
                   loc(G)==loc(Bwy) .or. loc(G)==loc(Bwz) .or. &
                   loc(G)==loc(ux) )then
                    G(0,:) = 0.0d0
                    do i=-1-ncx,-1
                        
                        G(i,:) = -G(-i,:)
                    end do
                else

                    do i=-1-ncx,-1
                        G(i,:) = G(-i,:)
                        
                    end do
                end if
            end if
        end if

        if(xMaxBoundary==opt)then
            if(ppi==ppx-1)then
                if(loc(G)==loc(ux) )then
                    G(Nx,:) = 0.0d0
                    do i=1,ncx+1
                        
                        G(Nx+i,:) = -G(Nx-i,:)
                    end do
                elseif(loc(G)==loc(Ex) .or. &
                    loc(G)==loc(By) .or. loc(G)==loc(Bz) .or. &
                    loc(G)==loc(Bfy) .or. loc(G)==loc(Bfz) .or. &
                    loc(G)==loc(Bwy) .or. loc(G)==loc(Bwz))then
                     
                     do i=1,ncx+1
                         
                         G(Nx+i,:) = -G(Nx-i+1,:)
                     end do
                else

                    do i=1,ncx+1
                        G(Nx+i,:) = G(Nx-i,:)
                        
                    end do
                end if
            end if
        end if

        if(yMinBoundary==opt)then
            if(ppj==0)then
                if(loc(G)==loc(Ey) .or. &
                   loc(G)==loc(Bx) .or. loc(G)==loc(Bz) .or. &
                   loc(G)==loc(Bfx) .or. loc(G)==loc(Bfz) .or. &
                   loc(G)==loc(Bwx) .or. loc(G)==loc(Bwz) .or. &
                   loc(G)==loc(uy) )then
                    G(:,0) = 0.0d0
                    do i=-1-ncy,-1
                        G(:,i) = -G(:,-i)
                    end do
                else

                    do i=-1-ncy,-1
                        G(:,i) = G(:,-i)
                        
                    end do
                end if
            end if
        end if

        if(yMaxBoundary==opt)then
            if(ppj==ppy-1)then
                if(loc(G)==loc(uy) )then
                    G(:,Ny) = 0.0d0
                    do i=1,1+ncy
                        G(:,Ny+i) = -G(:,Ny-i)
                    end do
                elseif(loc(G)==loc(Ey) .or. &
                    loc(G)==loc(Bx) .or. loc(G)==loc(Bz) .or. &
                    loc(G)==loc(Bfx) .or. loc(G)==loc(Bfz) .or. &
                    loc(G)==loc(Bwx) .or. loc(G)==loc(Bwz))then
                     
                     do i=1,1+ncy
                         G(:,Ny+i) = -G(:,Ny-i+1)
                     end do
                else

                    do i=1,1+ncy
                        G(:,Ny+i) = G(:,Ny-i)
                        
                    end do
                end if
            end if
        end if
    end subroutine useXinverseForGrid_in_exchange

    subroutine usePeriodicForGrid_Injector_for(G)
        implicit none
        real(kind=8) :: Tempx(2+ncx,-1-ncy:Ny+ncy+1), Tempy(-1-ncx:Nx+ncx+1,2+ncy)
        real(kind=8) :: G(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1)
        
        if(ppi==0 .and. xMinInjector)then
            if(yMinBoundary==1 .and. ppj==0)then
                G(0,0:Ny+ncy+1) = G(0,0:Ny+ncy+1)+G(inx,0:Ny+ncy+1)/2
            elseif(yMaxBoundary==1 .and. ppj==ppy-1)then
                G(0,-1-ncy:Ny) = G(0,-1-ncy:Ny)+G(inx,-1-ncy:Ny)/2
            else
                G(0,:) = G(0,:)+G(inx,:)/2
            end if

            G(-1-ncx:-1,:) = G(inx-ncx-1:inx-1,:)

        else if(ppi==ppx-1 .and. xMaxInjector)then
            if(yMinBoundary==1 .and. ppj==0)then
                G(Nx,0:Ny+ncy+1) = G(Nx,0:Ny+ncy+1)+G(Nx-inx,0:Ny+ncy+1)/2
            elseif(yMaxBoundary==1 .and. ppj==ppy-1)then
                G(Nx,-1-ncy:Ny) = G(Nx,-1-ncy:Ny)+G(Nx-inx,-1-ncy:Ny)/2
            else
                G(Nx,:) = G(Nx,:)+G(Nx-inx,:)/2
            end if
            
            G(Nx+1:Nx+ncx+1,:) = G(Nx-inx+1:Nx-inx+1+ncx,:)
            
        end if
        
        if(ppj==0 .and. yMinInjector)then
            if(xMinBoundary==1 .and. ppi==0)then
                G(0:Nx+ncx+1,0) = G(0:Nx+ncx+1,0)+G(0:Nx+ncx+1,iny)/2
            elseif(yMaxBoundary==1 .and. ppj==ppy-1)then
                G(-1-ncx:Nx,0) = G(-1-ncx:Nx,0)+G(-1-ncx:Nx,iny)/2
            else
                G(:,0) = G(:,0)+G(:,iny)/2
            end if

            G(:,-1-ncy:-1) = G(:,iny-ncy-1:iny-1)

        else if(ppj==ppy-1 .and. yMaxInjector)then
            
            if(xMinBoundary==1 .and. ppi==0)then
                G(0:Nx+ncx+1,Ny) = G(0:Nx+ncx+1,Ny)+G(0:Nx+ncx+1,Ny-iny)/2
            elseif(xMaxBoundary==1 .and. ppi==ppx-1)then
                G(-1-ncx:Nx,Ny) = G(-1-ncx:Nx,Ny)+G(-1-ncx:Nx,Ny-iny)/2
            else
                G(:,Ny) = G(:,Ny)+G(:,Ny-iny)/2
            end if

            G(:,Ny+1:Ny+ncy+1) = G(:,Ny-iny+1:ncy+Ny-iny+1)
        end if
        
    end subroutine usePeriodicForGrid_Injector_for

    subroutine usePeriodicForGrid_Injector_in_exchange(G)
        implicit none
        real(kind=8) :: G(-1-ncx:Nx+ncx+1,-1-ncy:Ny+ncy+1) 
        
        if(ppi==0 .and. xMinInjector)then
            if(loc(G)==loc(By) .or. loc(G)==loc(Bz) .or. &
            loc(G)==loc(Bfy) .or. loc(G)==loc(Bfz) .or. &
            loc(G)==loc(Bwy) .or. loc(G)==loc(Bwz) .or. loc(G)==loc(Ex))then
                    G(-1-ncx:-1,:) = G(inx-ncx:inx,:)
                    
            else
                    G(0,:) = (G(0,:)+G(inx,:))/2
                    G(-1-ncx:-1,:) = G(inx-ncx-1:inx-1,:)
            end if
        else if(ppi==ppx-1 .and. xMaxInjector)then
            if(loc(G)==loc(By) .or. loc(G)==loc(Bz) .or. &
            loc(G)==loc(Bfy) .or. loc(G)==loc(Bfz) .or. &
            loc(G)==loc(Bwy) .or. loc(G)==loc(Bwz) .or. loc(G)==loc(Ex))then
                    G(Nx+1:Nx+ncx+1,:) = G(Nx-inx+1:Nx-inx+ncx+1,:)
            else
                    G(Nx,:) = (G(Nx,:)+G(Nx-inx,:))/2
                    G(Nx+1:Nx+ncx+1,:) = G(Nx-inx+1:Nx-inx+1+ncx,:)
            end if  
        end if
        
        if(ppj==0 .and. yMinInjector)then
            if(loc(G)==loc(Bx) .or. loc(G)==loc(Bz) .or. &
            loc(G)==loc(Bfx) .or. loc(G)==loc(Bfz) .or. &
            loc(G)==loc(Bwx) .or. loc(G)==loc(Bwz) .or. loc(G)==loc(Ey))then
                    G(:,-1-ncy:-1) = G(:,iny-ncy:iny)
            else
                    G(:,0) = (G(:,0)+G(:,iny))/2
                    G(:,-1-ncy:0) = G(:,iny-ncy-1:iny)
            end if
        else if(ppj==ppy-1 .and. yMaxInjector)then
            if(loc(G)==loc(Bx) .or. loc(G)==loc(Bz) .or. &
            loc(G)==loc(Bfx) .or. loc(G)==loc(Bfz) .or. &
            loc(G)==loc(Bwx) .or. loc(G)==loc(Bwz) .or. loc(G)==loc(Ey))then
                    G(:,Ny+1:Ny+ncy+1) = G(:,Ny-iny+1:Ny-iny+1+ncy)
            else
                    G(:,Ny) = (G(:,Ny)+G(:,Ny-iny))/2
                    G(:,Ny+1:Ny+ncy+1) = G(:,Ny-iny+1:Ny-iny+1+ncy)
            end if
        end if
    end subroutine usePeriodicForGrid_Injector_in_exchange

end module Boundary
